// multinomial expression divergence
// (c) Yongjin Park, 2013
#include "util.hh"
#include "med.hh"
#ifdef MEDGAUSS
#include "med_gauss.hh"
#endif
#include <string>
#include <boost/lexical_cast.hpp>

template<typename T>
T get_param(const char* arg, std::string emsg)
{
  T val;
  try
  {
    val = boost::lexical_cast<T>(arg);
  } catch (boost::bad_lexical_cast&)
  {
    err_msg(emsg);
    exit(1);
  }
  return val;
}

////////////////////////////////////////////////////////////////
struct prog_args_t
{
  explicit prog_args_t() :
    data_file(""), ncols(0),  // required args
    tol(1e-3), max_iter(100), // options
    bmin(-2.), bmax(2.)
  {
  }

  std::string data_file;
  // std::string output;
  std::vector<size_t> ncols;
  double tol;
  int max_iter;
  double bmin;
  double bmax;
};

void print_help(const char* prog, const prog_args_t& args)
{
  std::cerr << prog << " -d data.txt -c c1 -c c2 ..." << std::endl;
  std::cerr << "         -bmin bmin -bmax bmax" << std::endl;
  std::cerr << "         -i iter -t tol" << std::endl;
  std::cerr << std::endl;
  std::cerr << " data.txt  expression matrix (n x p tab-separated)" << std::endl;
  std::cerr << " c[N]      # columns of N-th condition" << std::endl;
  std::cerr << "           such that sum_i c[i] = p" << std::endl;
  // std::cerr << " out       output files' header (default: " << args.output << ")" << std::endl;
  std::cerr << " bmin   min value of beta" << std::endl;
  std::cerr << " bmax   max value of beta" << std::endl;  
  std::cerr << std::endl;
}

bool parse_args(const int argc, const char* argv[], prog_args_t& args_out)
{
  for (int j = 1; j < argc; ++j)
  {
    std::string curr = argv[j];

    if (curr == "-d" && ++j < argc)
    {
      args_out.data_file = argv[j];
    }
    // else if (curr == "-o" && ++j < argc)
    // {
    //   args_out.output = argv[j];
    // }
    else if (curr == "-c" && ++j < argc)
    {
      args_out.ncols.push_back( get_param<size_t>(argv[j], "parsing ncol") );
    }
    else if (curr == "-i" && ++j < argc)
    {
      args_out.max_iter = get_param<int>(argv[j], "parsing iter");
    }
    else if (curr == "-t" && ++j < argc)
    {
      args_out.tol = get_param<double>(argv[j], "parsing tol");
    }
    else if (curr == "-bmin" && ++j < argc)
    {
      args_out.bmin = get_param<double>(argv[j], "parsing beta min");
    }
    else if (curr == "-bmax" && ++j < argc)
    {
      args_out.bmax = get_param<double>(argv[j], "parsing beta max");
    }

  }

  if (args_out.tol <= 0)
    return false;
  if (args_out.data_file.size() == 0 || args_out.ncols.size() == 0 )      
    return false;
  if (args_out.bmin >= args_out.bmax)
    return false;
  return true;
}


////////////////////////////////////////////////////////////////

#ifdef MEDGAUSS
typedef func_expr_pairwise_div_gauss_t func_type;
#else
typedef func_expr_pairwise_div_t func_type;
#endif

int
main(const int argc, const char* argv[]) 
{
  prog_args_t args;

  if( ! parse_args( argc, argv, args ) )
  {
    print_help( argv[0], args );
    exit(1);
  }

  const char* data_file = args.data_file.c_str();
  std::vector<size_t>& ncols = args.ncols;

  expr_tuple_vec_t data_vec;
  if( ! read_expr_tuple( data_file, ncols, data_vec ) )
  {
    TLOG("cannot read data");
    exit(1);
  }

  expr_tuple_set_t data_set;
  for(size_t j=0; j<data_vec.size(); ++j)
    data_set.insert( data_vec.at(j), 1. );

  func_type F(ncols.size(), args.bmin, args.bmax);

  const size_t max_iter = args.max_iter;
  const double tol = args.tol;

  for(size_t iter=0; iter<max_iter; ++iter)
  {
    // double t = (double) (iter + 1);
    F.update( data_set );
    double delta = empirical_bayes( F );
    TLOG( "iter = " << iter << ", penalty = " << F.get_penalty() << ", delta = " << delta );
    if( delta < tol ) break;
  }

  F.write_beta( std::cout );

  return 0;
}
